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ABSTRACT 

On a daily basis, consensus theory attracts more and more researches from 
different areas of interest, to apply its techniques to solve technical problems in a way 
that is faster, more reliable, and even more precise than ever before. A power system 
network is one of those fields that consensus theory employs extensively. The use of the 
consensus algorithm to solve the Economic Dispatch and Load Restoration Problems is a 
good example. Instead of a conventional central controller, some researchers have 
explored an algorithm to solve the above mentioned problems, in a distribution manner, 
using the consensus algorithm, which is based on calculation methods, i.e., non 
estimation methods, for updating the information consensus matrix. 

Starting from this point of solving these types of problems mentioned, 
specifically, in a distribution fashion, using the consensus algorithm, we have 
implemented a new advanced consensus algorithm. It is based on the adaptive estimation 
techniques, such as the Gradient Algorithm and the Recursive Least Square Algorithm, to 
solve the same problems. This advanced work was tested on different case studies that 
had formerly been explored, as seen in references 5, 7, and 18. Three and five generators, 
or agents, with different topologies, correspond to the Economic Dispatch Problem and 
the IEEE 16-Bus power system corresponds to the Load Restoration Problem. 
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In all the cases we have studied, the results met our expectations with extreme 
accuracy, and completely matched the results of the previous researchers. There is little 
question that this research proves the capability and dependability of using the consensus 
algorithm, based on the estimation methods as the Gradient Algorithm and the Recursive 
Least Square Algorithm to solve such power problems. 

The form and content of this abstract are approved. I recommend its publication. 

Approved: Miloje Radenkovic 
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1. Introduction 

For more than half of a century, the existing US power grid has been providing its 
electric power to a wide range of customers. In fact, this system faced substantial 
magnitudes of challenges that were not taken into account at the time the design was 
originated. These challenges can be summarized in accordance with the following: 

• lack of supporting the distributed generation. 

• renewable energy sources which are being added to the power grid. 

• lack of flexibility and adaptability. 

• natural disasters. 

• human operating errors. 

For all of the above reasons, and to assure that the power grid will continue to 
provide the demanded power to its consumers as needed, the power grid should be 
efficient and smart enough to deal with all demeanors of ongoing and varying changes of 
the system topology without losing its stability, availability, et cetera. 

In order to make a power system network effective and smart enough, it should 
implement a communication network and a control network that integrate and 
complement the future power system, where the system is called the Large-Scale 
Network Control System (LSNCS) [4] [5]. There is no doubt that the LSNCS will be a 
necessity for the next series of power networks. This refers to a specific means and 
existence of a large number of subsystems. In this bulk system, the conventional control 
schematic will face severe challenges that will inevitably affect its controllability. These 
challenges can be summarized in the following factors [5]: 

• requirement of a high level of connectivity to achieve the optimal operation 
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condition. 

• sensitivity to failures and modeling errors. 

• affect of a variety of the configurations, such as Plug-and-Play of the power grid 
topology which will make the smart grid topology unknown or undefined. 

To overcome the drawbacks of conventional control, distribution control is more 
convenient for solving the identified problems. The main task for using distribution 
control algorithms is to make sure all system agents, or nodes, reach what is called the 
consensus. Several areas, where the consensus algorithm is applied, can be used in the 
improvement of the power system. For example, the consensus algorithm of a multi- 
vehicle system, where a vehicle is able to communicate with other vehicles upon a 
communication network system, is a good example [2]. Therefore, to control this 
LSNCS smart grid, a very strong algorithm has to be implemented in order to work 
correctly in the existence of a lack of or unreliable communication capability; this applies 
even with loss of the central control topology. Such an algorithm has to be tested in areas 
that resemble LSNCS, as would the Economic Dispatch Problem and the Load 
Restoration Problem. 

To illustrate the idea of the central and distribution control, let us consider the 
"five buses power system", where each bus has its load and generated power as it is 
shown in Figure 1.1 on the next page. 
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(a) (b) 
Figure 1.1: (a) Central controller, and (b) Distributed controller 

The communication topology could be a conventional central controller, as shown 
in Figure 1.1(a), or a distributed controller as in Figure 1.1(b). In Figure 1.1(b), generator 
one has been selected as the consensus leader. The local controller, in each bus, will keep 
reporting its status to the leader. Based on the collected information, the leader will 
decide whether to increase or decrease the group of the consensus variable, according to 
the situation of a required condition [5]. The entire process mentioned is manipulated by 
the consensus algorithm. 

The consensus algorithm is used to solve some fundamental problems of the 
power system, as the Economic Dispatch problem and the Load Restoration problem. 
Instead of the conventional central controller, Some researchers had explored consensus 
algorithms to solve the mentioned problems in a distribution manner [5] [8] [18]. The main 
reason of using the distribution fashion is to overcome the drawbacks of the central 
controller. 

Staring from the point of solving the economic dispatch, and the load restoration 

problems in a distribution fashion, by using the consensus algorithm, we have 

implemented a new advanced consensus algorithm based on the adaptive estimation 

techniques of the Gradient Algorithm and the Recursive Least Square Algorithm to solve 

the same problems. This advanced work is tested on different case studies that are 
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proposed in the references [5] [7] [18] to ascertain the sets of comparatively- sourced 
results do coincide and are depicted accurately. These case studies, and their results, are 
presented in details in the Chapter 4. No doubt, this research proves the capability and 
dependability of using the consensus algorithm, based on the estimation methods, citing 
specifically, the Gradient Algorithm and the Recursive Least Square Algorithm to solve 
such power problems. 

Briefly, our research objectives outlining this thesis are summarized in the 
following section: 
1.1 Thesis Objectives 

To fulfill the goals of this research, a set of objectives was specified inclusive in this 
work. These objectives are summarized in the following points: 

1. Implementing a Matlab code to solve the Economic Dispatch Problem for the 
previous case studies that were focused on in the reference [5] without the power 
constraints and fulfilling the same results. 

2. Implementing a Matlab code to solve the Economic Dispatch Problem and taking 
into account the power constraints in order to fulfill the power limitations, i.e., a 
condition that was mentioned in Reference [7]. 

3. Implementing a new Matlab code to solve the Economic Dispatch Problem in the 
same case studies addressing objectives 1 and 2; but this time by applying a new 
estimation method using the Gradient Algorithm and fulfilling the same results in 
objectives 1 and 2. 

4. Implementing a new Matlab code to solve the Economic Dispatch Problem in the 
same case studies in objectives 1 and 2; but this time by applying a new 



estimation method using the Recursive Least Square Algorithm and fulfilling the 
same results in objectives 1 and 2. 

5. Repeating objectives 3 and 4 to study the case where the system topology changes 
with time. 

6. Implementing a Matlab code to solve the Economic Dispatch Problem to study 
the case if some agents are found to be disconnected. 

7. Implementing a Matlab code to solve the Load Restoration Problem for the case 
study that was presented in reference [18] and fulfilling the same results. 

8. Implementing a new Matlab code to solve the Load Restoration Problem in the 
same case in objective 7; but this time by applying the new estimation method 
using the Gradient Algorithm and obtaining the same results in objective 7. 

1.2 Thesis Organization 

The remainder of this thesis is organized and described in methodologies as 
follows: 

• Chapter 2: provides a very important background about the graph and consensus 
theories in the networks; the consensus mathematical model is presented in this 
chapter 

• Chapter 3: presents one of the very useful applications of the consensus theory, 
which is solving the economic dispatch and the load restoration problems in a 
distributed manner; the mathematical model of the new consensus estimation 
method is presented in this chapter, as well 

• Chapter 4: presents the study cases, and their respective results, of the estimation 
methods to solve the economic dispatch and the load restoration problems; also, it 
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contains notes and comments in direct reference and concern to the results 
• Chapter 5: thesis conclusions and recommendations for future studies are 
presented in this chapter, which is the last chapter of the thesis 
In the following chapter, the thesis starts with an introduction about the graph and 
the consensus theory. 
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2. Graph and Consensus Theory 

Before going in the details of the consensus theory, there is an important concept 
which must be clearly understandable. This concept is called the Graph. 
2.1 Graph Theory 

If we look around us in the real world, there are many situations that can be easily 
described or converted to a diagram which consists of a set of points, called vertices, 
together with lines, called edges, which connect or join certain pairs of the points within 
the diagrams. For instant, the subway system, where the points could represent the stop 
stations, and the edges represent the rails that connect these stations. The mathematical 
concept of a situation of such a case provides the idea of what is called a Graph. 

The principle of the Graph was born from the idea of the Konigsberg Bridge 
Problem. The Pregel River in Prussia divided the city of Konigsberg into four separate 
regions which were connected by seven bridges as shown in Figure 2.1 below. 




Figure 2.1: Bridges of Konigsberg city (1736) 
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A simple diagram of the bridges in the city is illustrated in Figure 2.2 below. 
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Figure 2.2: The Four regions and seven bridges in the Konigsberg city 
The inhabitants of the city wondered if it were possible to leave home, cross each 
of the seven bridges only once, and return home. The Leonhard Euler, a Swiss 
mathematician, (1707-1783), took the puzzle and tried to apply a mathematical method to 
solve the problem. He redrew each region as a vertex and each bridge as an edge 
connected vertices corresponding to the regions, as shown in Figure 2.3. 




Figure 2.3: Eulerian Circuit 
Figure 2.3 shows the Graph that encodes all necessary information from Figure 
2.2. It is called the Eulerian circuit. Many consider the Euler' s method to be the birth of 
the Graph theory [1][9][11]. 
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2.1.1 Graph Definition 

The Graph G is defined as a structure which consists a vertex set and an 

edge set £"(G), and a relation that associates with each edge whose two vertices are called 
the endpoints. This graph is used to model the system topology of the network. The graph 
G = (V, E) where V = {1,2, ...,nj is called the vertices or the nodes. E V x V is called 
the edges which is a set of unordered pairs of distinct vertices. When the two vertices 1, 2 
in V(G) are endpoints of an edge, both 1 and 2 are called adjacent [1]. 

In the real world, Graph G can be an electric circuit where vertices represent 
sources, diodes, capacitors, et cetera, and edges represent the wires that connect them. It 
can be a computer network where vertices represent computers, and edges represent the 
communication network that connects them together. It can be the World Wide Web, 
where vertices represent the WebPages, and edges represent the hyperlinks, and so on. In 
this research, the graph is an electric power grid where the vertices represent the power 
generators or loads, and the edges represent the communication network that connects 
them as illustrated in Figure 2.4 below, specifically, for the three generators case study. It 
is an unordered, an undirected, and un-weighted graph. 




Figure 2.4: Three generators (agents) system 
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2.1.2 Directed and Undirected Graph 

As was mentioned, Graph G is a mathematical structure consists of a set of 
vertices V(G), and a set of edges E(G), which connect the vertices according to an 
associated relationship. This mathematical structure can be a directed or an undirected 
graph. In the undirected graph, the edges that connect the vertices are not directionally 
fixed. The direction form a node n t to a node tij is represented by a line, with or without, 
arrowhead in both line ends, as seen in Figure 2.4 on the previous page. On the other 
hand, in the directed graph, also referred to as the digraph, the edge that connects two 
vertices is directed, or ordered, to gravitate in a specific direction from a node n„ or the 
tail node, to a node tij, called the head node. The edge direction can be easily recognized 
by the arrowhead that always points to the head node as illustrated in Figure 2.5 below. 




Figure 2.5: Directed graph 
Formally; G = (V,E) where V = {1,2, ...,n} and E QV XV . 
Assuming that u and w are end points then: 

G — (V E) will be ( unc ^ recte ^ if f or allu,w E V , i. e (u, w) ££"<=> (w, u) E E 
^directed otherwise 



10 



2.1.3 A Connected and Disconnected Graph 

The Graph G is a connected graph if, between every two vertices in the V(G), 
there is a path linking them, otherwise the Graph G is called a disconnected graph. In 
other words, the graph is disconnected if its vertex set can be divided into two non-empty 
subsets where no edge connects them [10]. Figure 2.6 below illustrates the two types. 




(a) (b) 
Figure 2.6: (a) A connected graph, and (b) a disconnected graph 



2.1.4 An Incidence and Adjacency Matrices 

In the graph G, when the edge e joins or connects the two vertices 1 and 2, these 
two vertices are called the adjacent. In the same time, both 1 and 2 are called the incident 
to the edge e. Likewise, in the graph G, two edges are adjacent if they have a vertex in 
common as it is shown in Figure 2.7 below. 




(a) (b) 
Figure 2.7: (a) Adjacent vertices, and (b) adjacent edges 
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Thus, based on the previous description, a pertinent question is raised. Are the 
graphs, as they are in diagram form, suitable for the computer calculations or the 
mathematical models to study their properties? The answer is no. Because of this, these 
two matrices are considered critical matrices. They are the incidence matrix and the 
adjacency matrix. Both of them are associated with the graph. 

• The incidence matrix of a graph G is n x m matrix M G := (m ve ), where m ve is the 
number of times (0,1, or 2) that a vertex v and an edge e are incident. 

• The adjacent matrix of a graph G is n x n matrix A G := (a uv ), where a uv is the 
number of edges that join vertices u and v. 

Both the incidence and the adjacency matrices of the graph G in Figure 2.6 (a) can be 
built as follows: 



It is easily apparent to notice that the adjacency matrix is much smaller than the 
incidence matrix because usually graphs have more edges than vertices. This causes the 
adjacency matrix to require less storage computer capability which is preferred in the 
computer programs where the consensus algorithm is applied [10]. 



The incidence matrix; M- 



7x11 — 



1 000110000 0" 

11000010000 

01100000001 

00110000110 

00011001000 

00000101100 

000001001 1- 



The adjacency matrix; A 7x7 = 



oiooiio- 

1 1 1 
10 10 1 
10 111 
10 10 10 
10 110 
1110 
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2.2 Consensus Theory 

As mentioned, the concept of the consensus is an area of interest from many 
researchers in different fields of a study. To illustrate the idea of the consensus, let us 
consider a system of three decision-making units where each unit can be known as an 
agent, or a generator as shown in Figure 2.4. Each agent has two parts: first part is the 
local controller and the second part is the consensus manager. The local controller 
collects system data and reports its own condition, or state, to the consensus manager. 
The mission of the consensus manager is to analyze this information and negotiates it 
with other connected neighborhood agents via the communication system. As result, the 
consensus manager will pass the decision back to the local controller. The local controller 
will update its situation depending on the received information and report back its new 
condition to the consensus manager again in a reciprocal, or back-and-forth process. 

2.2.1 Consensus Definition 

In a network of multi-agents, a "Consensus" means to reach an agreement 
regarding a certain quantity of an interest, depending on the state of all agents. A 
"Consensus Algorithm", or a protocol, is an interaction rule that specifies the 
information exchange between an agent and all its neighbors in the same network [2]. 

2.2.2 Consensus in Networks 

The topology of the system in Figure 2.4 can be represented using undirected 
graph G = (V,E) where V = {1,2, ...,n} and E Q V x V, whereas this graph can be 
described as using the adjacency matrix A as mentioned prior, of a finite system of n 
nodes. The adjacency matrix A is n x n matrix that depends on the number of the agents. 
The off-diagonal entry a,y is the number of edges from node i to node j. If a y - ± , i j 
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then the agent i th is communicating with the agent f h . In our case of a simple finite graph, 
the adjacency matrix is only a (0, l)-matrix with zeros on its diagonal. The neighbors of 
i th agent are denoted by: JVJ = {j G V ■ (i,j) EE). 

For our system, the i th node has X; G R that denotes the state value it possess. The 
state Xj of a node can be a physical quantity, such as voltage, current, power, incremental 
cost, et cetera. It can be said that the system nodes, generators, or agents, have reached 
the consensus state if, and only if, the state value of i th node is equal the state value of j th 
node for all i , j [3]. i.e 

Xj = Xj for all i ,j (2.1) 
The information discovery process for agent i is presented as follows: 

Xi k+1 = x t k +I ljE K t a i] (xj k - Xi k ) , i = 1,2,... n (2.2) 

Where: 

Xi k+1 ,Xi k : are the discovered information by agent i (like average net power) 
j\f; : is the neighbors that are connected to the agent i 

a t j : is the coefficient for information exchange between neighboring agents i and j 
Oij = I 

The information discovery process for the whole system is formed using matrix format as 
follows: 

X k+i =x k + A*X k (2.3) 
X k+1 = (/ + A)X k (2.4) 
X k+1 = DX k (2.5) 



_ (0 < a t j < 1 if agent i and j are connected 
otherwise 
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Where: 



yk _ r v k v k v k~\ 
L 1 ' ■■■ ' i ' ■■■ ' 71 J 



D = 



1 ^}£N 1 



a. 



a 



tl 



a 



it 



a 



lit 



i] 



a 



The D matrix is called the information discovery matrix, or sparse iteration matrix, for 
the system [18]. 

The coefficients a i; in matrix D have to be properly designed, because they play a 
critical role in reaching consensus of a system. For instant, a fully-connected system can 
reach consensus in one step if a i; is selected to be equal - where n is the number of the 

system nodes. However, most of the systems are usually not fully-connected where they 
can suffer from loss of their connection lines. 

Different methods are used to compute matrix D based on computing the 
coefficients a i; . Some of those methods are summarized in the following: 

• The Uniform Method: the weights are fixed and computed according the equation 
(2.6). The drawbacks of this method are the loss of the adaptivity and the slow 
convergence. 



f i 

n 



CLa — < 



y otherwise 



(2.6) 



• The Metropolis Method: it is an improved method to accommodate the changes of 
system configuration which guarantees the stability of the information discovery 
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process. The weights are updated according the equation (2.7) below. 

1 



[max (rii ,nf)+ 1] 



l_y.^ I i = j (2-7) 

L l eK i [max (n ( ,^)+l] 7 

^0 otherwise 
• The Mean Metropolis Method: a new improved method that is presented in [18]. 
It provides the stability and the adaptivity for the information discovery process. 
Their weights are updated according the equation (2.8) below. 

j G N t 



[n ; + rij + s] 

^0 otherwise 



1 - V 2 / - / (2.8) 



Where: 

s: is a very small number, and can be zero for large complex systems. 
2.3 Consensus Algorithm Applications 

Here are some applications where the consensus algorithm plays an important 

role:- 

• Flocking Theory: A group of moving agents that have sensing and 
communicating devices in a large environment. The mission of the consensus 
algorithm in here is to make every agent reach the same velocity of its neighbors. 
The system of the network is a dynamic system with a dynamic topology. The 
topology is an agent state dependence, and it is computed locally for each agent. 
A good example for this is a group of non-pilot planes [2] [12]. 

• Rendezvous in Space: In this case, the position of the agent is very important 

which makes the system an interactive position topology. This application is 

easier than the Flocking Theory because it does not require an agent-to-obstacle 
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clash avoidance [13]. 

Fast Consensus in Small -Worlds: Using semi-define convex programming, where 
the weight of the network is considered, leads to a slight increase in algebraic 
connectivity of the network, where the speed of the convergence of the consensus 
algorithms is measured. By keeping the network weights fixed, the topology can 
be designed to achieve a relatively high algebraic connectivity. 
Distributed Sensor Fusion in Sensor Networks: The above applications are 
considered as distributed sensor fusion. Distributed averaging problems require 
implementing the Kalman filter. The average of the inputs in the sensor networks 
is dynamically computed [14]. 

Distribution Formation Control: In multivehicle systems, the design and analysis 
framework of distributed formation controllers uses vectors of relative positions 
of neighboring vehicles. In order the agents, as multivehicle, to move in 
formation, the agents have to cooperate, which requires consent and collaboration 
of every agent in the formation. It is clearly a consensus problem [16]. 
Economic Dispatch and Load Restoration Problems in Power System: a 
conventional centralized control problem can be solved in a distribution manner 
by applying the consensus theory. Surly, this requires to choose the appropriate 
measures as a consensus variables [5] [8] [18]. More details are presented in the 
next chapter. 
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3. Economic Dispatch and Load Restoration Problems 
3.1 Economic Dispatch Problem 

The main target here is to get perfect use of serving equipment as much as 
possible with parallel to the highest financial benefits. In a power system, the Economic 
Dispatch is an optimal operation point of each power generator where the demand power 
is met, but the total generation cost is the lowest [6]. 

Several methods are used to solve the Incremental Cost equations, like setting 
value of Lambda and solving for generation values, or solving the value of Lambda 
directly for a specific load demand with different algorithms. Modern techniques use 
linear and non-linear programming like Piecewise-Linear fuel-Cost functions and 
Quadratic programming of the cost function. 

Usually, a cost fuel curve of any power generator is formulated by quadratic 
function as follows [7]: 
The fuel cost equation: 

Ci(P Gi )= a i +p i P Gi + YiPgi (3-D 
The Total fuel cost equation: 

Crotai = (3-2) 

Under the balance condition: 

Pd ~ XUPa = (3.3) 

Where: 

Cj : the fuel cost of the i th generator 
a i>Pi>Yi '■ the cost coefficients of the i th generator 
P D : the total demanded power 
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P Gi : the generated power of the i generator 

The definition of the Incremental Cost IC for each generator is the same as the 

conventional economic dispatch: 

Ki = dc^) = x _ . = 12 n (34) 

By using the discrete-time consensus algorithm form then the Incremental Cost by using 
the first order discrete consensus algorithm is as follows [5]: 

h [k + l]= Z"=i dy Aj [k] , i = 1,2, ... , n (3.5) 

Where: 

dij : is the ( i ,j ) entry of the row-stochastic matrix D n depends on the Laplacian 
matrix 

A t [k + 1] : updating IC of i th takes in account the connected j th generator to the i th 
generator 

In order to include the power balance constraints in the equation (3.3), the power balance 
error AP is defined as follows: 

&P= P D - Z? =1 P Ci (3.6) 
Finally, the updating IC of the leader becomes as follows: 

h[k + 1] = £ ; n =i d i} Aj [k] + sAP , i = 1,2, ...,n (3.7) 

Where: 

e : is a positive scalar, and it is called the convergence coefficient that controls 
the convergence speed of the leader generator or the leader agent. 
In the case where the power generation constraints are applied, the following conditions 
must be taken in account in the algorithm [5] [7] [8]. 
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Pd — ^Gi max ,whenP Gi > Panax 

Pa = Pa > when Pa_min ^ Pa ^ Pa max 
Pa = Pa mm .when Pa < Pa min 

3.2 Load Restoration Problem 

One of the consequences, after the protection system gives its orders to the circuit 
breakers to open, in order to isolate a faulty area in the power network, some un-faulted 
loads may trip and become out of service. Returning these loads back to the service as 
fast as possible, is one of the important aspects in the power networks. This action takes 
into account priority and capacity of each load. This procedure (retuning un-faulted loads 
back to the service) is called in the power system the "Load Restoration" [18] [19]. 

This kind of control task has been achieved in past years by using different 
computing algorithms, like genetic algorithm, particle swarm optimization, et cetera. A 
major drawback of these methods appears to be their need to a powerful central controller 
in order deal with a vast amount of information throughout the network. This leads to an 
extra cost and suffers from a single-point failure besides other drawbacks of a central 
controller that were mentioned earlier. In absence of question, solving the load restoration 
problem in a distribution control scheme is an intelligent solution to overcome all central 
controllers' drawbacks. 

Recently, Multi-Agent System (MAS) is an area of the interest, and different 
algorithms are applied to solve the load restoration problem. Some of them work with 
only radial power structures; others work with ring power structure. One advanced 
algorithm is presented in the reference [18]. It is capable to work even with a complex 
power structure, i.e., radial, ring, or mixed, because the proposed algorithm is 
independent of a system configuration. 
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In this case, the total net power of each agent is required for the calculation, and 
during stable normal operation the average of the total net power must be equal to zero. 

X eq = S^ii = o, i = l,2...n (3.8) 

Where: 

X : is discovered information which is the total net power of the agent i th 

n : is the total number of the agents 
The overall information updating process i.e., the total average net power for each agent, 
is modeled the same as the equation (3.7) but without the balance error dimension. 

X[k + 1] = DX[k] (3.9) 

Where: 

D : is the sparse iteration matrix. 

During a post fault period, the total net power equation (3.8) cannot be fulfilled 
because of the negative power shortage from the un-faulted agents, or loads. This 
requires immediate load restoration procedure. The Optimal load restoration decision can 
be obtained by satisfying coming two conditions: 

1. Maximizing the capacity of restored loads and assuring higher priority loads be 
served first, i.e., 

maxY.UPi.Pu (3.10) 

Where: 

Pj : is the priority index associate with i th load 

P Li : is the amount of power for an agent i th needs to be restored 

2. Sum of the total loads to be restored must not exceed the total generated power, 
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i.e., 




Ill Pu ^0 



(3.11) 



Where: 



P Gi : is the amount of generated power for an agent i . 

Identifying lost loads to be restored can be easily done by using the 0-1 Knapsack 
Problem, which is a combinatorial optimization problem. Having a group of items where 
each item has its own weight and corresponding value, and finding the number of each 
item to be added in a collection with a condition of the total collection weight, is less than 
or equal to a given limited amount but the total value is as large as possible. Its 
mathematical model is summarized as follows [21]: 



Pi : is the value of an agent i 1 

Wj : is the weight of an agent i th 

A : is the limited total weight of the collection 

The 0-1 knapsack problem can identify restored loads based on three factors, 
greedy by profit, or value, greedy by weight, or greedy by profit density. Since these 
loads are identified, the consensus algorithm will be activated to achieve a power balance 
again where the total net power for each system's agent should converge to zero. 
3.3 Estimation Information Matrix in the Consensus Algorithm 

Instead of computing the information matrix D in equations (3.5) and (3.9) based 
on the calculation methods that are presented in the chapter two, the estimation adaptive 



r 



maxf(x 1 ,x 2 ,...,x n ) = YH=\Vi x i 
^ =l w iXi <A 
x t £{0,1} 



(3.12) 



v 



Where: 
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techniques are applied to estimate it. In this thesis, Matlab codes were implemented based 
on the Gradient Algorithm (GA) and the Recursive Least Square Algorithm (RLS) to 
solve both the economic dispatch and the load restoration problems. The final forms of 
those estimation methods are presented in the following sections. The Complete 
mathematical derivation for both estimation algorithms is presented in the appendix A. 
3.3.1 Gradient Algorithm Estimation Form 

A t (k + 1) = diMk) + eN . dn [Aj(k) - Ai(k)] (3.13) 

du=l- SyeMjdi; . i =1,2,..., n (3.14) 
Inspired by the adaptive control theory, we propose the following cost function to be 
minimized [22]. 

min^i^C/c + l) - Xtik + 1)} 2 (3.15) 

dij (k + 1) = dij (fc) - ii[Ai<ik + 1) - Uk + 1)] (3.16) 

d ij (k + l)= di.jik) - n[Ai(k + l) - Xi(k + l)][Aj(k) - Ai(k)] (3.17) 
d ij (k + l)= dij(k) - n<Pi(k)[Ai(k+l) - Xt(k + 1)] (3.18) 

<Pi(k) = [AjQi) - Uk)] j = 1,2 Nt (3.19) 

In the equation (3.18), X { is the average of the neighbors of a i th node, and it is computed 
as follows: 

Uk + 1) = fl ;EMi l ; (/c + l) (3.20) 

The initial values: d i; (0) = i j and (0) = 1 — e M . d i; (0) i j 

IX = 0.5 , 0.1 or less is a scalar to increase the speed of the convergence. 
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3.3.2 Recursive Least Square Algorithm Estimation Form 

dtj (fc + 1) = dij (k) - Pi(fc)^(fe)[Ai(fc + 1) - Mfc + 1)] (3.21) 
<^(/c) = [XjQi) - Ai(fc)] 7 = 1,2 Nj (3.22) 

In the equation (3.21), 1; is the average of the neighbors of a i th node, and it is computed 
as follows: 

Aj(k + l) = ^-Z ;eN ^-(fc + l) (3.24) 

The initial values: d i; (0) = i 7 and 0^(0) = 1 — 2; e N t dij (0) i 7 

Pi(0) = Pj/ , Pj > is a scalar matrix to increase the speed of the convergence. 
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4. Study Cases and Results 

To verify our new estimation methods of computing average consensus, in both 
cases, the Economic Dispatch and the Load Restoration problems, they were tested on 
cases that had been studied before. We compared our results with theirs to make sure the 
performance of our new estimation algorithms gave the same correct results of the 
previous studies. 

4.1 Three Generators System for Economic Dispatch Problem 

The system contains three generators, or agents, which provide power supply to 
an electric load P D . The system parameters are summarized in the table 4.1 below [7]: 



Table 4.1; Three generators; system parameters 



Agent 
Or 
Generator 


Cost Coefficients 
used in equation (3.1) 


Power initial 


a 
[$/hr] 


P 

[$/MWhr] 


[$/MW 2 hr] 


[MW] 


1 st 


561 


7.92 


0.001562 


450 


2 nd 


310 


7.85 


0.001940 


300 


3 rd 


078 


7.97 


0.004820 


100 



In this case, the new estimation method is applied to solve the economic dispatch 
problem. Here, instead of computing the row- stochastic matrix D, based on the 
Laplacian Matrix L, the gradient algorithm is used to estimate it. The communication 
topology of this system is the same as Figure 1.1 (b). The algorithm will solve the 
economic dispatch problem so that the Incremental Cost of each generator converges to 
the same optimal value. The next sections show our results based on the estimation 
method. 
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4.1.1 Three Generators System with New Estimation Method using Gradient 
Algorithm without Power Constraints 

In this case, there were no produced power limits for the system generators. The 
situation's summary is as follows: 

1- The system parameters were taken from the table 4.1 

2- The row- stochastic matrix D was estimated by using the Gradient Algorithm 

3- The required demanded power was 850 MW and 950 MW 

4- The convergence gradient gain coefficient \i — 0.5 
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(a) (b) (c) 

Figure 4.1: Three generators using GA; IC's and power (no power constraints) 

As it is see in Figure 4.1 (a), all the three Incremental Costs X x ,X 2 , and X 3 

converged very fast to the same optimal value X* = 9.1483 %/MWh. After 1 second, the 

demanded power increased to 950 MW. We notice that the leader (agent 1), red line, 

responded faster than the others after the change happened, then it raised the other agents 

to a new Incremental Cost value, and after 1.4 sec they reached to the consensus again, 

which is a new optimal value X* — 9.295 %/MWh. In both stages, the generated power 

reached the demanded power quickly. Because there were no power constraints, we 

notice that all agents produced power without any limits as it is shown in Figure 4.1 (c). 
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Also the power balance had been satisfied as it is shown in Figure 4.1 (b). Our results are 

exactly the same as the results in the reference [5] where the Incremental Cost was 

computed based on the calculation methods in equations (2.6), (2.7), and (2.8). 

4.1.2 Three Generators System with New Estimation Method using Gradient 
Algorithm with Power Constraints 

In this case, there were produced power limits for the system generators. The 
situation's summary is as follows: 

1. The system parameters were taken from the table 4.1 except = 0.00128 
$/MW2hr so that agent one can produce more than 600 MW 

2. The row- stochastic matrix D was estimated by using the Gradient Algorithm 

3. The required demanded power was 850 MW and 950 MW 

4. The convergence gradient gain coefficient \i — 0.5 

5. The generator maximum's limits were [ 600 400 200 ] MW 

6. The generator minimum's limits were [ 150 100 050] MW 
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(a) (b) (c) 

Figure 4.2 Three generators using GA; IC's and power (with power constraints) 

In this case, the generators' power limits were applied in solving the Incremental 

Fuel Cost. As we see in Figure 4.2 (a), all the three Incremental Costs X 1 ,X 2 , and X 3 
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converged very fast to the same optimal value X* = 8.576 $/MWh. After 1 second, the 
demanded power increased to 950 MW. We noticed that the leader (agent 1), red line, 
responded faster than the others after the change happened, then it raised the other agents 
to a new Incremental Cost value, and after 1.6 sec they reached to the consensus again, 
which is a new optimal value A* = 8.852 $/MWh. In both stages, the generated power 
reached the demanded power quickly. Because of the power constraints, we notice that 
the generator one (the leader) stuck on 600 MW as a result of its power limits as shown in 
Figure 4.2 (c). Also, the power balance and the power constraints had been satisfied as 
shown in Figure 4.2 (b). Our results are exactly the same as the results in the reference 
[7] where the Incremental Cost was computed based on the calculation methods. 
4.2 Five Generators System for Economic Dispatch Problem 

The system contains five generators or agents provide power supply to an electric 
load P D . The system parameters are summarized in the table 4.2 below [7]: 



Table 4.2: Five generators; system parameters 



Agent 
or 

Generator 


Cost Coefficients 
used in equation (3.1) 


Power initial 
[MW] 


a 
[$/hr] 


P 

[$/MWhr] 


Y 2 
[$/MW 2 hr] 


1 st 


561 


7.92 


0.001562 


200 


2 nd 


310 


7.85 


0.001940 


250 


3 rd 


078 


7.97 


0.004820 


100 


4 th 


561 


7.92 


0.001562 


200 


5 th 


078 


7.80 


0.004820 


100 
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In this case, the gradient algorithm and the recursive least square algorithm are 
used to estimate the row- stochastic matrix D. The communication network of this system 
has different topologies as it is shown in Figure 4.3 below. The Incremental Cost 
algorithm will solve the economic dispatch problem so that the Incremental Cost of each 
generator converges to the same optimal value. The next sections show our results based 
on the estimation methods. 




Figure 4.3: Five generators; different system topologies 



4.2.1 Five Generators with New Estimation Method using Gradient Algorithm 

Three different cases for the initial values of the adjacency matrix for different 
system topologies were studied. The situation's summary is as follows:- 

1- The system parameters were taken from the table 4.2 

2- The row- stochastic matrix D was estimated by using the Gradient Algorithm 

3- The required demanded power was 850 MW and 950 MW 

4- The convergence gradient gain coefficient \i — 0.5 

5- Different system topologies were used as it is shown in Figure 4.3 

6- No power constraints were applied 
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4.2.1.1 Five Generators using GA; Guessing Initial Values dj: = djj = 0.0 



dij = dji = 0.0 , d a = 1 — ^ dj 



j=l:N t 



rl.O 


0.0 


0.0 


0.0 


O.Ol 


0.0 


1.0 


0.0 


0.0 


0.0 


0.0 


0.0 


1.0 


0.0 


0.0 


0.0 


0.0 


0.0 


1.0 


0.0 


LO.O 


0.0 


0.0 


0.0 


1.0-1 



Incremental Costs 



f 



0.5 



A - a 7565 











1 £- m 


561 







1 

Time [sec] 



1.5 



1000 



950 



900 



850 



Total Power (Demand & Produce! 



750 



700 



0.5 



1 1.5 
Time [sec] 



(a) (b) 
Figure 4.4: Five generators using GA; IC's and power (dy = djj = 0.0) 
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4.2.1.2 Five Generators using GA; Guessing Initial Values = djj 0. 



dij = dji ^ 0.0 , d u = 1 — ^ dy 
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Figure 4.5: Five generators using GA; IC's and power (d^ = djj 0.0) 
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4.2.1.3 Five Generators using GA; Guessing Initial Values djj ^= djj ^=0.0 



dij * dji * 0.0, 



d a = 1 



- X di 

j=l:N t ,t*j 
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Figure 4.6: Five generators using GA; IC's and power (dy djj 0.0) 

It can be seen in Figures 4.4 (a), 4.5 (a), and 4.6 (a) by applying the gradient 

algorithm, all five Incremental Costs X 1 ,X 2 A3 ,X 4 ,andX 5 converged very fast to the 

same optimal value X* = 8.666 $/MWh. After 1 second, the demanded power increased 

to 950 MW. It is noticeable that the leader (agent 1), red line, responded faster than the 

others after the change happened, then it raised the other agents to a new Incremental 

Cost value, and after 1.1 sec they reached to the consensus again, which is a new optimal 
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value A* = 8.756 $/MWh. It can be easily noticed that in Figure 4.5(a) where dy = dji ^ 
0.0, the agents reached the consensus faster than the other cases because of the symmetry 
of the initial matrix. Figures 4.4 (b), 4.5 (b), and 4.6 (b) show the total generated power 
by the agents. When the demanded power changed from 850 MW to 950 MW after 1 
second of operation, all agents increased their produced power to meet the demanded 
power. Therefore, the power balance had been satisfied. Our results are exactly the same 
as the results in the reference [5] where the Incremental Cost was computed based on the 
calculation methods. 

4.2.2 Five Generators System with New Estimation Method using Recursive Least 
Square Algorithm 

Three different cases for the initial values of the adjacency matrix for different 
system topologies were studied. The situation's summary is as follows:- 

1- The system parameters were taken from the table 4.2 

2- The row- stochastic matrix D was estimated by using the Gradient Algorithm 

3- The required demanded power was 850 MW and 950 MW 

4- Different system topologies were used as it is shown in Figure 4.3 

5- No power constraints were applied 

6- The Covariance Matrix's initial p = 10/ 
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4.2.2.1 Five Generators using RLS; Guessing Initial Values djj = dj t = 0.0 



dij = dji = 0.0 , d u = 1 — ^ dj 
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Figure 4.7: Five generators using RLS; IC's and power (djj = djj = 0.0) 









P D 

P G 






I 




[ 

























34 



4.2.2.2 Five Generators using RLS; Guessing Initial Values djj = d jt 0. 



d-ij = dji ^ 0.0 ,d u = l— ^ 
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Figure 4.8: Five generators using RLS; IC's and power (djj = dji 0.0) 
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4.2.2.3 Five Generators using RLS; Guessing Initial Values djj dj t 0. 
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Figure 4.9: Five generators using RLS; IC's and power (cL cL 0.0) 



It can be seen in Figures 4.7 (a), 4.8 (a), and 4.9 (a) by applying recursive least 
square algorithm, all five Incremental Costs Xi ,X 2 ,^3 , A, 4 , and A 5 converged very fast to 
the same optimal value A* = 8.666 $/MWh. After 1 second, the demanded power 
increased to 950 MW. It was noticeable that the leader (agent 1), yellow line, responded 
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faster than the others after the change happened, then it raised the other agents to a new 
Incremental Cost value, and after 1.05 seconds they reached the consensus again, which 
is a new optimal value X* = 8.756 $/MWh. It can be easily noticed that in Figure 4.8(a) 
where dy = dji 4- 0.0, the agents reached the consensus faster than the other cases because 
of the symmetry of the initial matrix. Overall, the performance of RLS algorithm has 
smoother convergence to the optimal value than the GA algorithm, and this can be easily 
noticeable when their results are compared. Figures 4.7 (b), 4.8 (b), and 4.9 (b) show the 
total produced power by all agents. When the demanded power changed from 850 MW to 
950 MW after 1 second of operation, all agents increased their produced power to meet 
the demanded power. So the power balance had been satisfied. Our results are exactly the 
same as the results in the reference [5] where the Incremental Cost was computed based 
on the calculation methods in equations (2.6), (2.7), and (2.8). 

4.2.3 Five Generators with New Estimation Method using GA and RLS algorithms 
where the System Topology Changes with Time 

In this case, the system topology changes with time. The situation's summary is 
as follows:- 

1 . The system parameters were taken from the table 4.2 

2. The row- stochastic matrix D was estimated by using the Gradient Algorithm and 
the Recursive Least Square Algorithm with do = 0.0 for both algorithms 

3. The required demanded power was 850 MW and 950 MW 

4. The system topology was changing with time as it is shown in Figure 4.10 on the 
next page 

5. No power constraints were applied 
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The topology of the system changes with time as follows:- 
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Figure 4.10: Five generators; changing sequence of the system topology 
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Figure 4.11: Five generators using GA; IC's and power (topology changes with time) 
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Figure 4.12: Five generators using RLS; IC's and power (topology changes with time) 

It can be seen from Figures 4.11(a) and 4.12(a), by using the gradient algorithm 

and the recursive least square algorithm, all five Incremental Costs X 1 , X 2 , A 3 , X, 4 , and A 5 

converged very fast to the same optimal value A* = 8.666 $/MWh and X* = 

8.756 $/MWh corresponding to the demanded power from 850 MW to 950 MW. Figures 

4.11(b) and 4.12 (b) illustrate the total generated power. It is noticeable that the total 

generated power increased to meet the demanded power. This case lead to the conclusion 

the connected agents reached the consensus even though their topology changed with 

time. 

4.2.4 Five Generators with New Estimation Method using GA and Agents 4 and 5 
are Become Disconnected 

In this case, two agents became disconnected from the rest of the system. What 

will happen to the system convergence and the produced power? The changing sequence 

of the system topology is illustrated in Figure 4.13 on the next page. The agents 4 and 5 

have been selected as the disconnected agents after 1.2 sec of operation. The situation's 

39 



summary is as follows:- 

1. The system parameters were taken from the table 4.2 

2. The row- stochastic matrix D was estimated by using the Gradient Algorithm 

3. The required demanded power was 850 MW 

4. The convergence gradient gain coefficient n — 0.5 

5. No power constraints were applied. 




(0.4 to 0,6") sec (0.6 to 1.2) sec (1.2 to 2) sec 

Figure 4.13: Five generators; agents 4 and 5 become disconnected 
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Figure 4.14: Five generators using GA; IC's and power (two agents 4 and 5 become 

disconnected) 

It can be seen from Figure 4.14 (a), by using the gradient algorithm, all five 
Incremental Costs X 1 ,X 2 ,^3 , ^4, and X 5 converged very fast to the same optimal value 
X* = 8.666 $/MWh. After 1.2 sec of operation, two agents (4 & 5) became disconnected 
from the rest of the system while the demanded power still at 850 MW. It can be noticed 
that the other three agents (1, 2, and 3) converged to a new Incremental Cost value 
X* = 9.148 $/MWh. As the system lost the connection with the disconnected agents so it 
kept the last Incremental Costs of them just before they became disconnected. Because 
the system lost part of the generated power from the disconnected agents, as it is seen in 
Figure 4.14 (b), the other agents had to compensate the shortage in the generated power, 
and this what exactly happened in a few milliseconds as it illustrated in Figure 4.14 (c). 
The agents 1, 2, and 3 increased their produced power to match the demanded power of 
the system. Finally, the power balance had been satisfied. 
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4.2.5 Five Generators with New Estimation Method using RLS and Check the 
Necessity of the Reaching Consensus 

This case is to check the necessity that said "In order the agents to reach the 
consensus, the sum of all columns of the adjacency matrix d must be equal to one", so the 
summation of all rows (red lines) and columns (blue lines) were plotted in order to check 
whether this necessity hold or not. The situation's summary is as follows:- 

1. The system parameters were taken from the table 4.2 

2. The row- stochastic matrix D was estimated by using the RLS and d D = 0.0 

3. The required demanded power was 850 MW 

4. No power constraints were applied. 

The system topology was changing with time as it is shown in Figure 4.15 below. 




(0.4 to 0.6) sec (0.6 to 2} sec 

Figure 4.15: Five generators; changing sequence of system topology 
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Figure 4.16: Five generators using RLS; IC's and power (P Q = 0.01) 
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Figure 4.17: Five generators; adjacency rows and columns' summation (P = 0.01) 
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Figure 4.18: Five generators using RLS; IC's and power (P G = 10) 
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Figure 4.19: Five generators; adjacency rows and columns' summation (P = 10) 
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It can be seen this necessity was accomplished when Po = 0.01 in Figure 4.17, but 
it did not hold as P D increased as in Figure 4.19. It is noticeable the system reached the 
consensus even though the sum of columns of the row- stochastic matrix D did not equal 
one. From this ironic result, it can be concluded that this proposed necessity has not to be 
held true in order for the system to reach the consensus. Instead, the average of the 
column's summation of the row- stochastic matrix D must be equal to one. 
4.3 IEEE 16-Bus Test Power System for the Load Restoration Problem 

In this case, the system contains 16 buses compound of generation and load 
agents. The system configuration is show in Figure 4.20 below: 
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Figure 4.20: The 16-bus test power system (source ref. [18]) 
It is a standard IEEE test power system. The arrows resemble the power loads, the 
circles resemble the power generators, the red squares resemble the circuit breakers in the 
close position, and the green squares resemble the circuit breakers in the open position. In 
this case of study, the power loses of the transmission lines are neglected. 
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Data of the 16-bus system in Figure 4.20 on the previous page can be easily summarized 
in the table 4.3 below [18]: 



Table 4.3: 16-Bus power system; system data 



Agent 
or 

DUs 


Neighbors 


[MW] 


[MW] 


[MW] 


Priority 


1 th 


2,5 


40 





40 


- 


2 th 


1,3,5,10 


40 





40 


- 


3 th 


2,10 


50 





50 


- 


4 th 


10,11 





25 


-25 


P2 


5 th 


1,2,6 





35 


-35 


P2 


6 th 


5,7,9 





20 


-20 


PI 


yth 


6,8,9 





15 


-15 


P2 


8 th 


7,9 





25 


-25 


PI 


9 th 


6,7,8,11,15 





5 


-5 


PI 


10 th 


2,3,4 





30 


-30 


P2 


11 th 


4,9,12 


40 





40 




12 th 


11,13 





45 


-45 


P2 


13 th 


12,14 





10 


-10 


P2 


14 m 


13,16 





40 


-40 


P2 


15 m 


9 


30 





30 




16 m 


14 


50 





50 





In the table above; P Gi is the agent i local generation, P Li is the agent / local 
load, X° is the agent i th local net power, and Pi~2 is the agent i th load priority with Pj > 
P2. The assumption here according to [20] is a sudden fault occurs to the generator on the 
bus 1 1 (agent 1 1) which led the protection system to react immediately to isolate the fault 
by opening some circuit breakers. The system configuration of the post fault period is 
illustrated in Figure 4.21 on the next page. 
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Figure 4.21: The 16-bus system; post fault configuration (source ref. [18]) 



As a result the power protection reaction, the system lost eight un-faulted loads 
(out-of- service) which attached to the agents 4 and 6-12 (shaded area). No doubt, these 
loads should be restored, as fast as possible, based on the total net power available and 
their priority. So, the load restoration process has to be activated to return back the lost 
loads. For such case, very important information like the total net power, the agents' 
indexes, and the demanded power of the un-faulted loads are extremely required. 

To solve the load restoration problem for the situation above in a distribution 
fashion, the information discovery matrix M is very important. Each agent has to 
converge to the same M matrix. According to the reference [18], every agent in the 
system is initially created with the nx3 matrix. In this initial matrix only two nonzero 
elements can be found. Mu indicates the agent index (i th ) if it is working normally or (0) 
if it is not. indicates the out-of- service agent index (i th ) if it is ready for restoration or 
(0) if it is not. Mjj indicates the local net power for the working agent or the demanded 
power of the ready for restoration agent. Clearly in the Mijs matrix, the i th agent row 
could be [i Pu ] if it is working agent or [0 i Pu ] if it is ready for restoration agent. 
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Because every agent has to converge to the same information matrix M, the 
information discovery algorithm should be applied to each element of the initial matrices. 
Here, instead of apply information discovery algorithm based on the calculation methods, 
which are proposed in the reference [18], a new advanced information discovery 
algorithm which is based on the estimation methods like, the gradient algorithm, is 
applied. The final converged matrix M contains all necessary information required for the 
load restoration. The first column shows all normally working agents, the second column 
shows all ready for restoration agents, and finally, the third column shows the local net 
power for the agents whether they are working normally or ready for the restoration as it 
is shown in the table 4.4 below. 

Table 4.4: 16-Bus power system; final converged of information matrix-M 



Initials matrices for the agents 



Final converged matrix M 



1 st column 



2 nd column 



3 rd column 



1 O 40 





ro 

o 

u 

ro 





o o - 



-25 








o 
2 40 




o J 

o 



O 7 

o 



15 





o 
o 

13 




O O 




5 





O S 

o o 



-35 








-25 




ro o 

3 o 

o o 

o o 

O 6 





50 

o 



— 20 



O 


O 9 



10 

o 

o o 



30 
O 



-o 





o- 















-0 





0- 





O —10 




o o 

14 O 



o 

-40 



J 






16 50 



O O 
O O 

12 



o 

15 









o 
o 



-45 

o 












30 









1/16 





40/16 


2/16 





40/16 


3/16 





50/16 





4/16 


-25/16 


5/16 





-35/16 





6/16 


-20/16 





7/16 


-15/16 





8/16 


-25/16 





9/16 


-5/16 





10/16 


-30/16 








0/16 





12/16 


-45/16 


13/16 





-10/16 


14/16 





-40/16 


15/16 





30/16 


16/16 





50/16 
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Using the final converged matrix M as input to the 0-1 knapsack problem, the 
exact ready for the restoration agents with their associated power can be easily identified 
to reconnect them back to the system as fast as possible. At the end, the balanced system 
is achieved. The final arrangement of our example is illustrated in Figure 4.22 below for 
the 0-1 knapsack greedy by value. The agents with red circuits were not selected by the 0- 
1 lknapsack greedy by value for the restoration. 




□ disconnected 



Figure 4.22: The 16-bus system; after restore configuration (source ref. [18]) 
The operation sequence in the Matlab code of this case was as follows: 

1- Total Number of iterations N = 800. 

2- Fault occurred at N = 300 

3- Load restored at N = 500 

4- The convergence gradient gain coefficient \i — 0.0001 
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4.3.1 IEEE 16-Bus Test Power System with New Estimation Method using GA and 
0-1 knapsack Value Based 
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Figure 4.23: The 16-bus system; information matrix-M convergence 
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Figure 4.24: The 16-bus system; total net power convergence (knapsack value based) 
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Instead of using the consensus algorithm based on the calculation methods that 
were proposed in the reference [18], a new advanced estimation method using the 
gradient algorithm was applied, and the 16-bus system reached to the same consensus 
results as that presented in the reference [18]. For instant, the information discovery 
matrix-M was initiated with initial values corresponding to every agent as in the table 4.4, 
then it converged to the final converged matrix-M shown in the same table. Figure 4.23 
shows the agent convergence process for each matrix column and Figure 4.24 shows the 
system total average net power convergence in the three stages process; the pre fault, the 
post fault, and the after fault period. In each stage, the system converged to the net power 
initial values. When the system was working normally (1-300 iterations in Figure 4.24) 
all agents converged to zero which is clear if we take the average of the agents in the 
column five of the table 4.3 on the page 46. During the post fault period (300-500 
iterations in Figure 4.24) and because of the out-of- service agents, remains working 
agents converged to a positive value which is 7.81 MW. Those working agents can be 
easily identified using the first column of the final converged matrix-M ((40+40+50-35- 
10-40+30+50)/ 16 = 7.81). Since the average total net power is positive, there is a 
necessity to restore some, or all out-of-service agents, depending on the system capacity, 
and their priorities. 

Finding such agents (their indexes) for the load restoration is very easy by 
extracting the information from the second column of the final converged discovery 
matrix-M. Sometimes and because of power demands' shortage, not all out-of-service 
agents can be restored. To identify agents that can be restored, the 0-1 knapsack is 
applied to perform such selection based on the agents' priorities and capacities. After 
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reconnection out-of- service agents as it is shown in Figure 4.24 (500-800 iterations), the 

total net power converged back to zero because the 0-1 knapsack value based restored 

100% of the demanded power available (7.81x16 = 125 MW). The 0-1 knapsack selected 

agents was; agent 4(-25), agent 8(-25), agent 10(-30), and agent 12(-45). 

4.3.2 IEEE 16-Bus Test Power System with New Estimation Method using GA and 
0-1 knapsack Density Based 
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Figure 4.25: The 16-bus system; information matrix-M convergence 
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Figure 4.26: The 16-bus system; total net power convergence (knapsack density based) 
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This time, the 0-1 knapsack density based was applied for the load restoration. 
By applying the 0-1 knapsack density based, only 96% of the demand power available 
(7.81x16x96% = 124 MW) was restored. The 0-1 knapsack selected agents was; agent 
4(-25), agent 6(-20), agent 7(-15), agent 8(-25), agent 9(-5), and agent 10(-30). As it is 
seen in Figure 4.27 (500-800 iterations), the average total net power converged to a small 
positive number which is (5/16=0.3) as a result of 96% of the 0-1 knapsack density 
based. 

Both sections 4.3.1 and 4.3.2 prove the capability of using our new advanced 
consensus algorithm based on the estimation methods like the gradient algorithm for 
solving such power problem. Moreover, they prove that it is needless to apply traditional 
calculation methods like those which are presented in the references [5] and [18]. 
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5. Conclusion and Future Work 
5.1 Conclusion 

In this thesis, a new concept of solving the economic dispatch problem and the 
load restoration problems is presented. As a part of solving both problems in a 
distribution fashion, an important information matrix is required for updating system 
status so a system can reach what is called the consensus. This matrix is called the 
information discovery matrix, and it is computed based on specific calculation methods. 

Instead of using the existing calculation methods of computing the information 
discovery matrix, adaptive techniques, like the gradient algorithm and the recursive least 
square algorithm, are applied to estimate above mentioned matrix. This advanced work is 
extremely useful to save time consumed during the calculation process of reaching 
consensus especially in the case of the bulk network systems. Our results of reaching the 
consensus totally matched previous results of cases that had used calculation methods of 
reaching the consensus. No doubt, this work proves the importance of the estimation 
methods in the reaching consensus and their capability of solving power problems in a 
distribution manner which helps to overcome the drawbacks of the central controller. 

In this thesis, chapter two gives an important background about the graph and 
consensus theory. The chapter three gives the theoretical analysis of applying the 
consensus algorithm to solve the economic dispatch and the load restoration problems. 
The estimation mathematical model is presented in this chapter. Our case studies and 
their results are presented in the chapter four. The comments and notes on these results 
are presented here also. Finally, the chapter five gives the conclusion and the future work 
of this research. 
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5.2 Future Work 

Some important steps have to be taken as further investigation to complement this 
work of study in the consensus algorithm based on the estimation techniques. These 
future researches can be summarized in the following points: 

• Study the case of a leader selection in a case if the system is a weighted, an 
ordered, and a directed graph. 

• Study the case of taking into account the power losses in the transmission lines of 
the power network. 

• Study deeply the case of the disconnected agents whether the disconnection 
happen in the communication network only, the power lines only, or both of them. 

• Study the case of losing the system's leader. Dose the system able to select 
another leader or not? And how? 

• Study the case of clastericity when the system is divided to several connected 
groups and each group converges to its own consensus. 

All in all, all these mentioned future researches need more of an advanced code that 
can take into account many possibilities and assumptions in order get a generally flexible 
algorithm. Such an algorithm can be easily applied in the real power networks. 
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APPENDIX A 

MATHEMATICAL DERIVATION OF THE ESTIMATION ALGORITHMS (RLS 

and GA) 

Let us consider the following linear system model [22]: 
y(i + 1) + a-^yii) + ••• + a n y(i + 1 - n) = b-^uii) + b 2 u(i - 1) + ••• + b m u{i + 1 — 
m),i = 0,1,2, ... . (A.l) 
Where: 

u(i) : measurable system input signal 

y(i) : measurable system output signal 
Knowing that q" 1 is unit delay operator, i.e. q _1 y(0 = y(i — 1) then system in equation 
(A.l) becomes: 

(1 + a x q~ x + ••• + a n q- n )y(i + 1) = + b 2 q~ 1 + - + b m q- m+1 )u(i) (A.2) 

y(i + 1)=^ U (0 (A.3) 
So that: 

y4((7 -1 ) = 1 + a^ -1 + ••• + a n q _n (A.4) 

Biq- 1 ) = b t + b 2 q~ x + - + b m q- m+1 (A.5) 

The term m equation (A.3) is transfer operator of the system. 

Let us define system parameters and signal vectors as follows: 

0l = [-a t , -a 2 , ... , -a n ; u(i), ... , b m ] (A.6) 

0(j) T = [y{i),...,y{i + 1 -n); ^, ...,u(i + 1 - m)] (A.7) 
So the system is equation (A.l) can be written as: 

y(i + 1) = 0(i) r o (A.8) 
More general system model is: 
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y(i + 1) = 0(Q T 9 O + d(i) (A.9) 
Where represents the cumulative effects of disturbance, measurement errors, etc. 
Assuming that system parameters are unknown, an estimate of o is denoted by 0(i + 1) 
then the following identification model can be constructed: 

y(i + 1) = 0(i + l) r 0(i) (A.10) 
Estimating 0(i + 1) is determined so that in a certain sense model output y(i + 1) 
matches the system output y(i + 1). 

Let us choose 6(i + 1) so that it minimizes the sum of the square errors of the fit. 

V(e) = IU (y(/c + 1) - 0(k) T 6(i + I))' (All) 

A necessary condition for a local minimum is: 

dv(e) . . 

^^ = which gives: 

2 SU (ytf + 1) - 0(k) T 0(i + 1)) (-0(fc)) = (A.12) 
It follows that: 

(IU 0(/c)0(/c) r )0(i + 1) = Ei=o 0(fe)y(fe + 1) (A.13) 
Hence a least square estimate of 6 is any estimate 0(i + 1) that satisfies equation (A.13). 
There may be more than one minimize of v(e). 
Let define a new matrix: 

R(O = n= o 0(k)0(ky (A14) 

Which is invertible, and called covariance matrix, then there is an unique least square 
estimate: 

0(i + 1) = ie(i)- 1 ZU 0(fe)y(* + 1) (A.i5) 
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A.l Derivation of Recursive Least Square Algorithm Estimation (RLS) 

Considering the least square estimate in equation (A. 13), we wish to determine a 
recursive form of this estimate that suitable for easily updating 0(0 to 0(i + 1). 
Recalling equation (A. 14), we see that 

R{i) = R{i - 1) + 0(i)0(O T (A.16) 
Then from (A. 13) 

R(i)6{i + l) = n =o 0(k)y(k + l) (A17) 
i?(i)0(t + l)= m0(fc)y(*+i) + 0(i)y(i + i) (A18) 
R(i)9(i + 1) = i?(i - l)0(i) + 0(i)y(i + 1) from (A. 13) (A. 19) 

R(Q§(i + 1) = ( R(0 - 0(O0(O r ) 0(0 + 0(Oy(i + 1) from (A.16) (A.20) 
Thus: 

R(i)9(i + 1) = R(i)0(i) - 0(O0(O r 0(0 + 0(Oy(i + 1) (A.21) 
R(i)9(i + 1) = i?(O0(O + 0(0 [y(i + 1) - 0(O r 0(0 ] (A.22) 
Then: 

0(i + 1) = 0(0 + i?(O -1 0(O[y(i + 1) - 0(O r 0(0 ] (A.23) 
From stored values of 0(0 and R(i — 1), and new observations y(i + 1) and 0(0, by 
using equations (A.20 & A.23) we can calculate the new values 0(i + 1) and R(i). 
Equations (A.16 & A.23) are called the Recursive Least Square Estimate (RLS). 
To avoid inverting matrix i?(0 -1 m equation (A.23), we use what called the matrix 
invariance lemma: 

(A+BCD)" 1 = A" 1 - A^C 1 + DA^B) 1 DA" 1 (A.25) 
Here, all matrices inverses on the right-hand side are exist. 
Let define P({) A then from (A. 17) 
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P(i)" 1 = P(i - l)" 1 + 0(i)0(O r (A.26) 
By applying the rule of matrix invariance lemma, we set: 

A = P(i)" 1 ,B = 0(/c) , C = 1 , and D = 0(k) T then the equation (A.25) becomes 

P(i) = P(i - 1) + Pa - 1)0(fe)0( ° rp °- 1) ( A 27) 

Then we have obtained the following RLS estimate that requires no inversion to be 
performed at each step. 

0(i + 1) = 0(i) + P(i)0(O [ y(i + 1) - 0(O r 0(O ] (A.28) 

P(i) = P(i -D+ Pa - 1)0(fe)0( ° rp °- 1) (A 29) 
P(0) = P(0) r > 

A.2 Derivation of Gradient Algorithm Estimation (GA) 

Given an estimate 6, the squared prediction error is: 

e(0 2 = i(y(i)-0(i-l) r 0) 2 (A.30) 
Its gradient with respect to is: 

V e e(i) 2 = -W ~ l)[y(0 - 0(i - l) r 0] (A.31) 

The direction of the negative gradient 0(i — l)[y(i) — 0(i — 1)0] is the direction in 
which a small change in will lead to a decrease in e(i) 2 . Thus one motivated to 
consider the algorithm 

0(0 = 0(i - 1) + M0(i - 1) (y(0 - 0(i - l) r 0(i - 1)) (A.32) 
Where [i > is a gain to increase the speed of the convergence. 
Equation (A.32) is called the Gradient Algorithm Estimate (GA) [22]. 
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A.3 Basic Flowchart for the Estimation Algorithms (RLS and GA) 
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Set 

System parameters 
Initial conditions 
Number of iterations 



1 .Compute X based on initial 

conditions 
2. Guess initial values for dij 
matrix 



Consensus algorithm 



Compute power produced 
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Updating dij matrix based on 
GA or RLS estimation's 
algorithm 




Plot results 



End 



Figure A.l: Flowchart of solving EDP based on new estimation algorithms 
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